import pandas as pd
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

d = pd.read_csv("datos_polarizacion.csv")

angulo = d["angulo"]
intensidad_malus = d["intensidad_malus"]
intensidad_eliptica = d["intensidad_eliptica"]
intensidad_circular = d["intensidad_circular"]

theta = np.deg2rad(angulo)

# Ajuste no lineal de Malus
def malus(theta, I0, theta0):
    return I0 * np.cos(theta - theta0)**2

popt, pcov = curve_fit(malus, theta, intensidad_malus, p0=[0.43, 0.03])

I0_fit, theta0_fit = popt
errores = np.sqrt(np.diag(pcov))
print("Ajuste No Lineal:\n")
print(f"I0 = {I0_fit:.2f} +- {errores[0]:.2f} W")
print(f"Theta0 = {theta0_fit:.2f} +- {errores[1]:.2f} rad")

theta_fit = np.linspace(min(theta), max(theta), 200)
theta_fit_deg = np.rad2deg(theta_fit)

plt.figure()
#plt.scatter(angulo, intensidad_malus, label="Datos")
plt.errorbar(angulo, intensidad_malus,
             xerr=1, yerr=0.01,
             fmt='o', label="Datos", capsize=3)
plt.plot(theta_fit_deg, malus(theta_fit, *popt), 'r', label="Ajuste Malus")
plt.xlabel("Ángulo (º)")
plt.ylabel("Intensidad")
plt.legend()
plt.title("Ley de Malus")
plt.show()

# Ajuste lineal I0
x = np.cos(theta)**2
y = intensidad_malus

coef = np.polyfit(x, y, 1)
pendiente, ordenada = coef

print("-"*25)
print("Ajuste Lineal:\n")

coef, cov = np.polyfit(x, y, 1, cov=True)
pendiente, ordenada = coef

error_pendiente = np.sqrt(cov[0,0])

r = np.corrcoef(x, y)[0,1]

print(f"I0 = {pendiente:.2f} +- {error_pendiente:.2f} W")
print(f"r = {r:.5f}")

x_fit = np.linspace(min(x), max(x), 100)

plt.figure()
#plt.scatter(x, y, label="Datos")
plt.errorbar(x, y, yerr=0.01,
             fmt='o', label="Datos", capsize=3)
plt.plot(x_fit, pendiente*x_fit + ordenada, 'r', label="Ajuste lineal")
plt.xlabel("cos²(θ)")
plt.ylabel("Intensidad")
plt.legend()
plt.title("Ajuste lineal I0")
plt.show()

# Intensidad Elíptica
print("-"*25)
print("Intensidad Elíptica:\n")
plt.figure()
#plt.scatter(angulo, intensidad_eliptica)
plt.errorbar(angulo, intensidad_eliptica,
             xerr=1, yerr=0.01,
             fmt='o', capsize=3)
plt.xlabel("Ángulo (grados)")
plt.ylabel("Intensidad")
plt.title("Polarización elíptica")
plt.show()
e = ((min(intensidad_eliptica))/(max(intensidad_eliptica)))
print(f"excentricidad = {e:.2f}")

# Intensidad Circular

plt.figure()
#plt.scatter(angulo, intensidad_circular)
plt.errorbar(angulo, intensidad_circular,
             xerr=1, yerr=0.01,
             fmt='o', capsize=3)
plt.xlabel("Ángulo (grados)")
plt.ylabel("Intensidad")
plt.ylim(0.15, 0.42)
plt.title("Polarización circular")
plt.show()